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ABSTRACT 


This thesis examines the performance of three different orbital propagators to determine 
which provide the best performance for use in Low Earth Orbit Rendezvous. The performance 
evaluation is based upon the propagator’s accuracy and the amount of time required to produce a 
solution. A Cowell high-fidelity propagator is used as a base line for comparison with an Encke 
and Clohessy-Wiltshire propagator. To further enhance the examination a Jacchia-70 
atmospheric model and a GEM-9 Geopotential model are used to provide perturbing acceleration 
inputs to the propagators. All comparisons are performed in a Local Vertical, Local Horizontal 
Reference Frame with the target spacecraft at the coordinate center. Tainting of the input data set 
by a prior processor make the findings suspect. Findings support the prediction that while the 
Cowell propagator is the most accurate it also takes the most time to achieve results. Also, the 
Clohessy-Wiltshire, while taking the least time is the most inaccurate. The Encke propagator 
deliveries the most balanced result. 
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I. INTRODUCTION 


The intercept of two large objects in a gravity free environment can either be very 
easy or very hard depending on whether or not you want both objects intact after the 
meeting. If not, a high energy collision works every time. For rendezvous in space, the 
object is for two spacecraft to meet in space, work together and separate with no damage 
to either. This requires a great deal of knowledge and planning involving a medium in 
which human beings have little intuition. To counter this shortcoming tools have 
developed to help. 

Orbital propagators have been used for years to predict the positions and 
velocities of objects in orbit. They are an integral part of space navigation. Their only 
real shortcoming is their size and complexity. A good high fidelity propagator can 
involve thousands of lines of code and many minutes of processing time in order to 
predict a spacecraft’s orbital location. While this is acceptable for mission planning, it is 
not very useful for real-time propagation. 

In September of 1993 NASA launched STS-51 which experimented with a new 
navigational tool. The Detailed Technical Objective 700, consisted of a GRID laptop 
computer with software written by NASA and a group of students from the Naval 
Postgraduate School (Ref. 1). As a precursor to an automated navigational system, this 
experiment was designed to demonstrate how on-orbit orbital prediction and graphical 
displays could assist astronauts navigating of the Shuttle Orbiter. Such an aid was found 
to be especially valuable during Proximity Operations and Rendezvous (Ref. 1). The 
programs RPOP and NPSLite flew and were used during an actual rendezvous. 

The problem with RPOP or NPSLite as bases for an automated navigation 
program lays in their design. RPOP used a set of differential equations developed by 
Clohessy-Wiltshire for orbital prediction while NPSLite used a “f ’ and “g” function type 
propagator (Ref. 1). Neither of the propagators took into account perturbations caused by 
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external forces except those inherent in the Two Body Problem. This was not an arbitrary 
decision. Normally the equations used to determine orbital position for a two body 
problem are not that complex, and therefore fairly simple for a computer to handle. The 
addition of various gravitation and atmospheric models to more accurately reflect small 
perturbing accelerations can greatly increase both program size and processing time. In 
the case of RPOP and NPSLite, it was assumed that two objects in nearly identical orbits 
are perturbed relatively the same. Therefore, when dealing with a purely relative 
reference frame these perturbations could be ignored with little loss of accuracy. The 
purpose of this thesis is to determine the validity of that assumption. 

As a result of the success of RPOP and NPSLite, NASA has decided to develop a 
flight qualified version of the RPOP software. It is NASA’s intention to have the 
software available for use during the future Space Shuttle - MIR mission in 1995. One of 
the difficulties in the flight certification is verifying that the algorithm upon which the 
software is based is sufficiently accurate. By expanding this query, we can quickly find 
ourselves asking which of the many orbital propagation models and algorithms available 
is the best for the purpose at hand. This is the second objective of this thesis. 

As previously mentioned, the on-orbit computation time a propagator requires can 
be critical. In order to be useful for the flight crew, a propagator should accurately predict 
the Orbiter’s position and velocity at least one orbit, approximately 90 minutes, in 
advance. Because even the most accurate propagator is vulnerable to errors, this 90 
minute requirement must be broken up into a series of steps typically on the order of a 
few seconds. This in turn means that the propagator must provide thousands of 
predictions for each 90 minute propagation. Additionally, to be of any use on-orbit, the 
software needs to perform all of these predictions repeatably every few seconds. For a 
high fidelity, time intensive propagator, this time restriction may render the algorithm of 
little value. 


2 



Because of the time constraint, the best propagator is not necessarily the most 
useful. The best propagator balances accuracy with low computation time. An ideal 
propagator would be best in both categories, but as will be seen this is rarely if ever the 
case. In fact, the ideal propagator may eventually depend on far more than just accuracy 
and speed, but for the moment this study will be limited to just these two criteria. It is the 
intent of this thesis to provide the initial look and possibly a guideline to future research 
in the ultimate choice of a best orbit propagation model and algorithm. 
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II. REFERENCE FRAMES 


Which reference frame a body is being observed in is integral to defining and 
describing the forces interacting with it. By definition, the Newton’s laws of motion 
apply only in an inertial reference frame. In order to translate to a non-inertial reference 
frame and use the same laws, one must translate those laws into an appropriate form. 

Integral to the definition of a reference frame is the coordinate system used to 
define location within said frame. Basically, all a coordinate system does is determine the 
location of a point in space relative to another arbitrary point. While this may sound 
simple, it is extremely important because it allows for the comparison of two objects 
within the same reference frame. This comparison and the forces which act between and 
on the objects is normally the true goal of the investigation. 

Typically a reference frame and coordinate system are chosen in such a way as to 
simplify the problem. An example which is directly applicable to this thesis is the 
relative reference frame. In comparing two objects which are both in motion relative to 
an inertial reference frame it greatly simplifies the problem if a reference frame and 
coordinate system is fixed to one object. 

For this thesis three different reference frames/coordinate systems are applicable. 
Primarily among these is the inertial reference frame in which all physical laws are based. 
Since the two objects of interest are both orbiting the Earth it is natural to remove the 
rotation of the Earth about its axis as a simplification. This involves an Earth Centered 
Earth Fixed (ECEF) reference frame/coordinate system. The relative reference frame 
uses one of the spacecraft as the center of the frame in order to eliminate any identical 
orbital motion between the two crafts. In all three cases transformation have been 
developed in order to translate between any of the coordinate systems without loss of 
accuracy. 
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A. THE M-50 INERTIAL REFERENCE FRAME 


Inertial by definition is a reference frame which is fixed in space or is moving 
with a constant velocity. The Inertial Reference Frame used by this thesis was the NASA 
standard M-50 frame. The choice of this frame was obvious because the data supplied by 
NASA for this thesis was given in terms of these coordinates. It is displayed in Figure 
2.1, as a three dimensional, rectilinear system defined by x, y and z axes. The x axis is 
defined by the Earth’s Vernal Equinox or the First Point of Aries. The z axis coincides 
with the Earth’s rotation vector, and roughly points toward the star Polaris. The y axis is 
defined by the cross products of the x and z axis unit vectors. 


Polaris 



Figure 2.1 Inertial Reference Frame 

The M-50 reference frame is Earth centered but the Earth does rotate relative to 
the frame. Because the Earth also is in motion about the Sun it could be argued that this 
is not truly an inertial reference frame. While this is true, a Geocentric Inertial Frame is 
sufficiently accurate for this investigation as long as the speeds involved are not very 
large relative to the speed of light, because then that Newtonian mechanics must be 
replaced with Einsteinian. 

B. WGS-84 EARTH CENTERED EARTH FIXED REFERENCE FRAME 

In order to eliminate unnecessary complications caused by the Earth’s rotation an 
Earth Centered Earth Fixed reference frame is used. The use of this reference frame was 
also necessitated by the fact that the GEM-9 Gravity Model uses the WGS-84 coordinate 
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system. Again, as seen in Figure 2.2, the WGS-84 is a three dimensional, rectilinear 
system defined by x, y and z axes. The z axis is defined by the Earth’s spin axis while the 
x axis is in this case defined by 0° Latitude and Longitude. Again the y axis is defined by 
the cross product of the x and z unit vectors. The entire reference frame is in rotation 
about the z axis at a rate equal to the Earth’s rotation rate. 


North Pole 



Figure 2.2 Earth Centered, Earth Fixed Coordinate System 

C. RELATIVE MOTION REFERENCE FRAME 

The relative motion reference frame is where most of the work takes place. 
Because we are looking at two objects orbiting the Earth in near identical paths, much of 
the orbital motion can be ignored. It is logical to assume that forces acting on the two 
objects identically in an inertial or ECEF reference frame can also be ignored with little 
loss in accuracy when a relative system is adopted. Whether this assumption is true or 
not is the subject of this thesis. 

There are two potential coordinate systems. The first is called Local Vertical, 
Local Horizontal (LVLH) and is defined relative to the orbit of the spacecraft, while the 
second is called Roll, Pitch and Yaw (RPY), is defined relative to the spacecraft itself. 
Both have advantages or disadvantages depending on what is required. 

1. Local Vertical, Local Horizontal Coordinate System 

LVLH is again a rectilinear system with x, y, and z coordinates. The x axis is 
defined by the direction of motion while the z axis points radially outward from the Earth. 


7 


The y axis is defined by the cross product of the x and z axis. As seen in Figure 2.3 the 
horizontal plane is defined by the x and y axis and the vertical plane by the x and z axis. 
A LVLH coordinate system is in constant rotation relative to an ECEF coordinate system 
(typically about the y axis.) The magnitude of that rotation is equal to the spacecraft’s 
angular velocity about the Earth’s center. Since the LVLH coordinate system is based 
upon the object’s orbit and not it’s attitude and because the orbiting objects are small 
relative to the earth, they can normally be treated as point masses. 

y 



Figure 2.3 Local Vertical, Local Horizontal Coordinate System 

2. Roll, Pitch and Yaw Coordinate System 

The RPY coordinate system is identical to that of an aircraft’s with the x or roll 
axis pointing out the nose of the spacecraft, the z or yaw axis out the spacecraft’s belly 
and the y or pitch axis pointing out the starboard wing. See Figure 2.4 for a graphical 
representation of this coordinate system. This is a very advantageous system for any 
sensors or personnel located on the spacecraft, but since the spacecraft may be in rotation 
about its center of mass it is not ideal when discussing the motion between two orbiting 
objects. 
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Figure 2.4 Roll, Pitch and Yaw Coordinate System 


D. COORDINATE TRANSFORMATIONS 

As discussed previously, it is critically important to be able to translate 
coordinates and forces between the various reference frames in order it insure the 
accuracy of the Equations of Motion used. The translation process can be a complex one 
involving detailed math, resulting in the addition of many new forcing terms into the 
Equations of Motion. These pseudo-forces, while having no existence in inertial space 
are very important in a rotating reference frame. Fortunately, many of these terms cancel 
when dealing with relative motion. 

The fact that all coordinate systems used are rectilinear is convenient because 
transformation matrices between systems exist and are well known. Typically a 
transformation can be view as a series of three rotations about the coordinate axis such 
that one coordinate system is translated into the other. The rotation about a single axis 
can be seen in Figure 2.5 . 
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a l 

a3,b 3 

Figure 2.5 Axis Rotation 

Mathematically this rotation can be done by multiplying a vector by a rotational 
matrix to get the rotated coordinates. This rotational matrix is determined from the unit 
vectors of the two axis. 

“ A A A "] — 

a x -b x <V& 2 h cos0 sin (j) 0 

A C B = a 2 -b x a 2 -b 2 a 2 -b 3 = -sin0 cos 0 0 (2.1) 

AAA 

a 3 -b x a 3 -b 2 a 3 b 3 \_ 0 0 1 _ 

The rotational matrix in Equation 2.1 can then be multiplied with two other 
rotational matrices to get the final transformation. The order of the multiplication is 
important as are the order of the rotations, to insure that the x axis in the first coordinate 
system lines up with the x axis in the new coordinate system. 

'C E =[ I C A \ A C B \ B C E ] (2.2) 

Transformation matrices were generated and used in translations between all 
coordinate systems to insure proper accuracy of the data. 
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III. THE PROPAGATORS 


The purpose of any propagator is to accurately predict a spacecraft’s position and velocity 
at some future point in time. There are many different mathematical methods to produce these 
future vectors, each with their own advantages and shortcomings. While they may differ in 
application or procedure, all claim Newton’s Universal Law of Gravity, Equation 3.1, as their 
source. 


r, Gmim.2 

F = -— 

r 2 (3.1) 

where G is the Universal Gravitational Constant and r is the distance separating the point masses 
mi and m 2 . 

To properly apply this physical law to the case of orbital mechanics, some general 
assumptions must be made. The first is that the two objects of interest are point masses which is 
valid if we assume the bodies are either spherically symmetrical or relatively small. Secondly, it 
is assumed that these two point masses are the only ones in the system. In the case of an object 
orbiting the Earth, the assumption of two point masses gives a sufficiently accurate answer. 
Given these underlying assumptions and applying Newton's Third Law of Motion we obtain 
Equation 3.2 which is the equation of the force on an orbiting spacecraft. 


mr = — 


GMm _ 


(3.2) 


Dividing out the spacecraft mass we have Equation 3.3 which governs motion for the 
Restricted Two Body Problem. The restriction is due to the presupposed assumption that the 
Earth’s mass M » m which in the case of a spacecraft orbiting the Earth is definitely true. This 
assumption additionally fixes the coordinate system at the Earth’s center since the spacecraft will 
have very little effect of the system’s Center of Mass. 
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( 3 . 3 ) 


The zero to the right of the equal sign in Equation 3.3 results from our assumptions that 
there are only two bodies in the system providing accelerating forces. This equation, 
unfortunately, only approximately describes the motion of a spacecraft about the Earth. In reality 
there are small perturbations caused by other planets, atmospheric drag, the fact that the Earth is 
not really a point mass, etc. An expression of the form, 


" KJ1V1 _ 

r = - —r+d P (3.4) 

r 

where a p is the sum of all the perturbing accelerations, is more accurate. Whether to include 
certain accelerations in a p depends on the order of magnitude of these perturbing forces acting on 
the body, and on the time interval over which they are allowed to affect its motion. 

There are a plethora of possible solution methods used to solve the second order 
differential Equation 3.4 (Ref. 2). The distinguishing features between the different propagators 
is the manner in which this equation is solved. For the purposed of this thesis only three 
propagators are compared, they being the most popular. 

A. THE COWELL PROPAGATOR 

The Cowell propagator is the highest fidelity propagator used. It’s accuracy is due to the 
fact that it solves the complete unsimplifed differential equation. It begins by breaking the 
second order differentials into two coupled first order different equations. 

r = v 

i GM _ (3-5) 

v = -- r r+d p 

These equations are solved simultaneously using a Runge-Kutta 4th order integrator. A built in 
feature of the code decreases the time step used in the numerical integrator scheme when the 
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solution varies to rapidly over a short time interval. The step size is chosen in an effort to 
contain numerical truncation and rounding errors below a specified tolerance level. 

While the Cowell propagator is generally the most accurate it is also the slowest. This is 
because it must perform the numerical computations over a relatively large number of small time 
steps in order to maintain accuracy. In addition, the time required significantly increases when 
atmospheric and/or gravity models are used to increase the propagator’s accuracy. 

B. THE CLOHESSY-WILTSHIRE PROPAGATOR 

This propagator uses a set of linearized Equations of Motion for relative motion 
developed by Clohessy-Wiltshire. The derivation is taken from Reference 2 and Reference 3. 
The solution of the equations provides the position and velocity of a spacecraft, the chaser, with 
respects to another, the target. In order to define vectors between these two objects an 
orthogonol reference frame is attached to the target and moves with it. The coordinate system is 
a form of Local Vertical, Local Horizontal. 

The vector position of the chaser is given by 

r -r T +p (3.6) 

which is differentiated with respects to the M50 inertial coordinate system giving 

F = f T + p+2( (bxp)+(bxp+ebx( cbxp) (3.7) 

where (omega) is the angular velocity of the target about the center of the Earth. If we take 

F = g + d p (3.8) 

where g is the acceleration due to gravity and A is the acceleration due to external forces, 
Equations 3.7 and 3.8 can be resolved into x, y and z components. 
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x = -g— + Ax + 2coz + coz + 1 o 2 x 
r 


y - ~g—+Ay 
r 

-gf "— -j + A z + g T - 2(ox -qjx + co 2 z 


Assuming that the distance between the target and chaser is much smaller than the orbital radius 
of the target, i.e (p) 2 « r 2 x, the following approximations can be made 


r “T + ^J 

g~g T [i-— 

\ r T J 

x x 

~g~ * ~g T — 
r rr 

y y 

~g~~ ~ ~8t 
t rr 


Z + rr 


(3.10) 


which by substitution into Equation 3.9 yields 


x — —g —i- A x +2 ox + coz + cq 2 x 
rr 


y--g T -+ Ay 

rr 


(3.11) 


Z — 2 g T -1- A z + 2 (fix + (tix + cq 2 z 

rr 
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If we make the additional assumption that the target's orbit is circular, then 

( 0=0 
(0 = 

and Equation 3.11 becomes 

x= Ax + 2caz 

y = A y -( 0 2 y (3.12) 

z = A z -2qk+3cq 2 z 

Since y is uncoupled from x and z in Equation 3.12, it can be solved simply and separately. 
Assuming a solution of the form 

A 

y = Asintuf + 5 cos can —- 

of (3.13) 

A and B are determined from the initial conditions, giving the solution 

y 

y = y n cos ox +—oosinox h— f( 1-cosox) 

0) CO (3.13) 

A 

y = y o cos£Uf - y g ( (0$,m(iX )+—- sin cot 

co 

The equations for x and z are more difficult because they are coupled. We begin 
integrating the equation for x once and substitute the result into the equation for z. 

z + 2co(2coz + Axt + ci)-3co 2 z-A z = 0 (3.15) 

Assuming that the solution to Equation 3.15 is of the form 

z = A sin cot + B cos cot — ————(3.16) 

CO CO CO 
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it can be substituted into the equation for x dot resulting in: 


x = 2oA sin cot + 2coB cos cot - 3 A x t - 3 a + 2- 


(3-17) 


Integrating, 


x = -lAcosox + 2B sin cot — A x t 2 — 3 at + 2—t + c? 

2 co 


and applying the initial conditions we obtain 


(3.18) 


1 2 

A = — (i 0 i — Ax) 

CO CO 

r _ 2 _. 

Xt Xo ^ Za t 
CO CO' 

ci = Xo~2coi„ 

2 2 

ci — x,,^—(z»l — Ax) 
co co 


(3.19) 


Substituting and reordering terms, 


*=* +z [6((ii—sinfiafll+iT——3r +z 9 -(1-costur) +4 -|r 2 +4r(l-costar) +A —(tar-sintar) 
° L I co j 'L co J L 2 co' J L® 

r , . fsinarl A,, , 

v=yJcostar] +y - +-^-(l-coscur) 

l co j CO' 

z =x f—=-(l—costar)]+ zj4 -3 cos car]+ z„ +4 —^-(tar-sintar) +A -4(l-costar) 

° co ' L co J L co' J lor J 


:=ij4cos<xr-3]+^[6ffi(l-costaf)]+4[2sinar]+4. -^smaz-St + A -^-(1-costar) 


(3.20) 


v =y a cos cot— y (; (tasintat) -f—- sin tar 

CO 


, r . r 2, ,1 fsintar 

i: = i l-2sinmfl + c fStusinfflfl+z L cos6;ir ] + 4!— (1-costar) +A- 

" L J " l ' l CO J L co J 
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The Clohessy-Wiltshire equations are convenient since they give answers in LVLH 
coordinates. This implies, assuming the radius and velocity vectors are obtained by the chasing 
vehicle, that no coordinate transformation is required for orbit determination. A drawback the to 
the solution given above is that they are based upon the assumption that the target orbit is 
circular. This produces an inherent error in the algorithm whose magnitude varies with distance 
between the vehicles and eccentricity of the two orbits. The greater the distance and eccentricity 
the greater the magnitude. These disadvantages most likely make the Clohessy-Wiltshire 
Propagator the least accurate of the propagators used. 

C. THE ENCKE PROPAGATOR 

While the Cowell method is probably the most accurate propagation method it can be 
cumbersome and time consuming, especially if the perturbing acceleration is much smaller than 
the force of gravity. Because this is the case for an Earth orbiting spacecraft an alternative 
technique based upon integrating the vector difference between the Two Body Solution and the 
perturbed solution rather than the explicit state vector found in the Cowell Method. This is the 
idea behind the Encke Method. The particular version of the Encke method developed here is 
taken from Reference 4. 

The Encke Method assumes that the spacecraft is traveling in a conical path called 
osculating orbit. This osculating orbit is determined by the position and velocity vectors at a 
given time fusing the classical Two Body Problem. Because of the perturbing accelerations, the 
actual orbit diverges in time from the osculating orbit and these deviations are what is found by 
the method. These vector differences are then added to the osculating orbit position and velocity 
vectors to give the actual state vectors at a future time. 

As stated previously, the osculating orbit is derived using the restricted two body. The 
method employed in Reference 4 is based upon the so-called “f ’ and “g” functions. These are 
determined from the position and velocity vectors at a prescribed time t 0 , and can be used to 
generate the classical orbital element. The f and g functions are given in Equation 3 (Ref. 4). 
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/ = -[(cos(£-£ 0 )-l] + l 


^[sin(E-E o )-(E-E 0 )]+(t-O 


f = -^-sm(E- E 0 ) 

V 

£ = 1 - — [l~ cos(E - £„)] 


(3.21) 


Position and velocity vectors of the osculating orbit are given by: 


r*e = fr 0 + gv 0 
V OS c=fc+gVo 


(3.22) 


In the Encke formalism, these vectors are used in the differential equation for the 
difference vector (8): 


d 2 8 i il; /i ( r 3 V 

= 1_ ^r r + a P 

dt 2 rL rL \ r 3 J ' 


osc osc 


(3.23) 


To avoid difficulties in evaluating r, slight modifications to Equation 3.23 are made. 


Since 


^) = U*)+5(r) 


(3.24) 


we have 


1 _^ = _ F (^ = i_ ( i + ^ 


(3.25) 


where 


<5*(<5-2r) 


(3.26) 
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The function F(q) can then substituted into Equation 3.23. This revised second order differential 
equation can now be integrated to find 8 and its first derivative v. When these are substituted 
into Equation 3.27, the position and velocity vectors for the desired time are found: 


r = /L + 5 


v = v OJC + v 


(3.27) 


While this propagator is more accurate than the Clohessy-Wiltshire Propagator since it 
does not depend upon to circular orbits it does restrict itself, due to the method of determining 
the osculating orbit, to planar orbits. This means that greater error is likely in the out of plane 
direction, the y axis in the LVLH reference frame. Whether this error is significant is yet to be 
determined. 
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IV. THE PERTURBING ACCELERATIONS 


The perturbing accelerations which affect the motion of a spacecraft in Earth’s 
orbit are many. They are caused by physical phenomenon existing outside the orbital 
system, such as solar radiation, because as well as due to the fact that many of the 
assumptions made in developing the Restricted Two Body Problem are too simplistic. 
The simplifying assumptions inherent in the Two Body Problem, such as disregarding the 
existence of three body effects, cause errors in the predictions. Regardless of the source, 
these perturbing accelerations should be considered in the propagation techniques in order 
to insure accurate predictions. 

A. THE PERTURBATIONS 

When considering a Low Earth Orbit (LEO), just about every conceivable 
perturbing force has an effect. Each of them must be reviewed and either incorporated 
into the propagator or discarded depending upon the magnitude of the perturbation they 
produce. 

1. Solar Pressure. 

The perturbation known as solar pressure is caused by the impingement of solar 
photons on the surface of a spacecraft. One’s first impulse is to discard this as negligible. 
Before doing so, one should consider that in space there is no counter force, such as 
atmospheric drag, to resist such a force. Over a sufficiently long time interval this 
perturbing force could build up to a significant level. Near Earth, the Solar Pressure 
Constant is 4.7 x 10" 5 dynes/cm 2 . 

Typically Solar pressure effects the satellite’s orbital eccentricity, but as implied 
above, it typically takes time for it’s effects to be felt. The time frame over which the 
propagator must be accurate for the purpose of this research is one orbit, or approximately 
90 minutes, and therefore this effect can be reasonably neglected. 
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2. Magnetic 


The Earth’s magnetic field, like any magnetic field, can induce a force on a 
magnetized or a current carrying body. The magnitude of this force is dependent upon 
the size of the magnetic fields and/or current. Such a force is often used by spacecraft via 
magnetic torquers to actively introduce a change in the spacecraft attitude. Even if there 
were no active use of the magnetic field, a spacecraft, after a sufficient time, develops a 
slight magnetic field. 

The difficulties in using the Earth’s magnetic field for atitude control is really a 
boon for an orbital propagator because the constantly changing direction and magnitude 
which creates great difficulties in developing a lengthy, direction-stable force does not 
create a very coherent perturbing force. Additionally, because the Earth’s magnetic field 
is so weak, only a fairly large current in the orbiter creates an appreciable effect. For 
these reasons, magnetic perturbations are ignored in this thesis. 

3. Third Body Effects 

The reality of orbital mechanics is that there are several gravitating bodies rather 
that the two assumed in the Restricted Two Body Problem. The Moon, the Sun and the 
other planets all have a gravitational effect on an orbiting spacecraft. To accurately 
model these effects the system should be expanded from two bodies to three, four, etc. 
This was attempted for a three body system and was found to be analytically intractable 
except in very limited circumstances. Extensions to a system of four or greater bodies 
turns out to be of very limited value for practical matters involving orbiting satellites. 

One possibility is to treat these effects like any other perturbation. It has been 
found that the accumulated effects, while they do exist, are again insignificant over the 
short duration of a single orbit. Therefore these perturbations are also ignored. 
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4. Earth Gravity Harmonics 


The Restricted Two Body Problem assumes that the Earth is spherically 
symmetrical and therefore can be represented by a point mass. In reality it is a large 
semi-homogenous ball. For a point mass, it can be assumed that the force of gravity is 
constant at every point on a sphere encompassing that point mass. For a semi- 
homogenous ball, every point on that same sphere can experience a different gravitational 
force. In order to take this variation into account a model of this variation has been 
developed. In the case of gravity, the concept of a gravitational potential provides this 
model. 


By mapping the gravitational potential at all points around the Earth it is possible 
to determine the acceleration due to gravity at any such point. The acceleration in a given 
direction is simply the gradient of the potential in that particular direction. The external 
geopotential is based upon 


where: 




n 


K (sin <t>)(C nm cos mX + S nm sin mX) 


(4.1) 


|i = The Earth’s Gravitational constant 

r = The magnitude of the spacecraft radius vector 

a = The Semi-major axis of the Central ellipsoid (Earth) 

n,m = Degree and order of each term 

<\> = Geocentric latitude 

X = Geocentric longitude 

C nm , S nm = Gravitational coefficients 

P n m = Associated Legendre functions 

N = number of terms used 
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The Gravitational Coefficients can be broken into two categories, C„o , the zonal 
harmonics or the Cm, , the tesseral and sectorial harmonics. The zonal harmonics are 
also know as the J n coefficients of which J 2 is the most dominant. The harmonics define 
the shape and distribution of the Earth mass. 

The J 2 effect can have a significant effect on a spacecraft’s orbit in a short amount 
of time. For this reason a gravity model must be examined to see if it has any substantial 
effect. Also it must be determined to what fidelity, i.e. number of coefficients, the model 
should have. 

5. Atmospheric Drag 

The final perturbing force considered, is that due to atmospheric drag. Even as 
high as 1000 km the Earth has a measurable atmosphere. While it is true that it is 
essentially only an indistinct plasma, there are still particles which can impact a 
spacecraft and cause a change in velocity. Also, as the altitude decreases the atmosphere 
becomes more dense and therefore more of a factor. 

For LEO spacecraft, atmospheric drag can be the dominant perturbing force. The 
drag a spacecraft feels is a function of velocity, atmospheric density and the spacecraft’s 
ballistic coefficient. The velocity is constant for a circular orbit, but as discussed 
previously, orbits are never truly circular. For elliptical orbits velocity is constantly 
changing depending on the point along the ellipse. The fact that the atmosphere’s density 
is also variable as a function of both position and time further complicates the problem. 
The ballistic coefficient can be assumed constant since it depends upon the shape of the 
Orbiter which does not change. 

The fact that velocity and density are variable means that drag is constantly 
changing, typically in a very short time. As a result of this rapid change, Atmospheric 
Drag has the potential to significantly effect an orbital prediction. Velocity is easily 
gained from an osculating orbit model, but density often requires a very complex model. 
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B. THE MODELS 


For the purpose of this thesis only two perturbing forces are used, Earth 
Gravitational Harmonics and Atmospheric Drag. There are numerous models, each with 
their own advantages and disadvantages. Models for this thesis were primarily chosen 
because of their availability, but fortunately, both are also extremely accurate. 

1. The GEM-9 Geopotential Model 

All Geopotential Models are essentially the same. Equation 4.1 is used to generate 
the gravitational potential at a specific location. They differ in the values of the 
gravitational coefficients and the number used. For this thesis the GEM-9 model using 
coefficients from Table 3 of Reference 5 are used. The model has a maximum resolution 
of 30 x 30, and uses the WGS-84 coordinate system. This is an ECEF standard 
coordinate system. 

2. The Jacchia 71 Atmospheric Model 

There are a number of potential Atmospheric models. Some are purely 
exponential, while others are more complex, taking into account the date, time, latitude 
and longitude and the effect of the sun on the upper atmosphere. The Jacchia Models are 
of the latter category making them large but more accurate. This model was developed by 
Jacchia in 1970 and takes into account the 11-year and 27-year solar cycles. It has the 
added advantage of using the WGS-84 coordinate system. 
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V. COMPARING PROPAGATORS 


A. THE DATA SET 

The easiest way to judge the accuracy of any algorithm is to compare it to actual 
data from the phenomenon being modeled. If the algorithm results in a relatively good 
match with the data the model and algorithm are considered acceptable. The usual tack 
taken by a modeler, is to start with the simplest model available to describe a 
phenomenon, and include effects (or fine tune the model) until such a time as the model 
is acceptable. It may happen that further refinement is possible but perhaps unnecessary 
to produce a workable model and algorithm. If not, it is generally back to the drawing 
board for the individual who is developing the model and algorithm. 

The data set in this thesis comes from a rendezvous mission performed by the 
Space Shuttle, and consists of a series of state vectors for the Space Shuttle and a 
rendezvous target. Each state vector is comprised of the spacecraft’s orbital position and 
velocity as well as a time stamp when the measurement was made. The entire data set has 
3238 state vectors taken at approximately 3.5 seconds intervals. The data set, as seen in 
Figure 5.1, contains a record of two orbits during which time the Orbiter performs a 
rendezvous with the target vehicle. Unfortunately, it was determined near the end of the 
project that the data could have been processed using a high-fidelity propagator to fill in 
the gaps. This fact limits the possible accuracy to that provide by the high-fidelity 
propagator. It also could produce a more accurate prediction to appear less accurate. 

In Figure 5.1, the data is presented as x versus z, but since the x coordinate 
closure rate is constant, the figure also demonstrates orbital positions as a function of 
time. All state vectors are provided in M50 inertial coordinates. Additionally, because 
the data is of an actual rendezvous, there are numerous thmster firings performed by the 
Orbiter during the time interval covered by the data set. These firings produce anticipated 
errors in the propagators, but since they are essentially impulses they are also readily 
identifiable. 
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Figure 5.1 Data Set in LVLH Coordinates 



B. THE PROGRAMS 

Since the motivation for this thesis is the efficacy of various propagators in the 
development of the RPOP software, it seemed logical to retain the C++ coding of the 
actual candidates. The majority of the code used in this thesis is the NPSLite software 
developed by LT Lee Barker (Ref. 6), but numerous modifications and additional code 
has been written by the author. 

The Cowell propagator used was originally developed by Robert Gottlieb and 
Mike Fraietta of McDonnell Douglas Corporation for NASA. Some code modification 
were made by LT Lee Barker (Ref. 6) to make the algorithm work. Additional 
modifications were incorporated by the author into the C++ code written for this thesis, 
but no changes to the core algorithm were made. In order to properly compare the 
performance of the Cowell propagator with the Clohessy-Wiltshire propagator, both the 
orbiter and the target state vectors were predicted and converted to the LVLH reference 
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frame. This was to nullify any common errors made by the algorithms which might have 
occurred in both spacecraft predictions 

The Clohessy-Wiltshire propagator was originally written by Scott Tamblyn of 
NASA’s Johnson Space Center and Jack Brazzel of McDonnell Douglas. Modifications 
to this algorithm were made to incorporate perturbing accelerations. The perturbing 
accelerations include a Jacchia 70 atmospheric model and a GEM-9 gravity model coded 
for use by the Cowell propagator. By using the same perturbing accelerations as the 
Cowell, a potential source of error was eliminated. Both the orbiter state vector and the 
perturbing accelerations had to be converted from the inertial frame to the LVLH frame 
prior to use by the Clohessy-Wiltshire propagator. 

The Encke Propagator used the “f and g” functions developed by LT Les 
Makepeace (Ref. 1) and coded by LT Barker (Ref. 6) for NPSLite. A Runge-Kutta 4-5 
integrator which integrates Equation 3.23 was incorporated by the author. Gravitational 
and atmospheric models were once again taken from the Cowell propagator to maintain 
coherence in the project. As with the Cowell propagator, state vectors were predicted in 
inertial coordinates and then converted to LVLH for the purpose of comparison. 

Numerous difficulties were experienced in compiling and running the Encke 
propagator when the gravitational and atmospheric models were added. Without the 
additional models the Encke propagator ran well. Addition of the perturbing forces 
caused a divergence of the orbit predictions. Unfortunately there was insufficient time to 
track this error down, and therefore an interesting portion of the project could not be 
accomplished. 

C. PERFORMANCE 

The performance of the different propagators has not been judged solely on the 
accuracy of the predicted state vectors, but also on how that accuracy varied as the time 
step is increased. Additionally, the amount of computation time required to perform the 
propagation is compared. It was anticipated that in order to obtain the most accurate 
prediction, a high fidelity propagator would be needed to break up the prediction time 
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interval into a series of smaller time steps over which each propagator could be run. 
Ideally, the propagator requiring the shortest computation time over the largest acceptable 
time step would optimize performance. 

For the purpose of this project four time steps were chosen. The smallest step was 
based on the separation between consecutive state vectors in the data set, approximately 
3.84 seconds. In order to insure that actual state vectors would be available for 
comparison, the remaining three time steps were simply a factor two, five and ten times 
that base time step. This resulted in time steps of 7.68, 19.2 and 38.4 seconds. Average 
computation time was determined by dividing the total amount of time required to 
propagate all 3238 state vectors over all four time steps by the total number of state 
vectors propagated. 

D. FINDINGS 

The Cowell and Encke propagators proved to be the most accurate as seen by 
Figure 5.2 (a) and (b). The Encke Coordinate and Velocity Differentials were virtually 
identical to those of the Cowell method. Differentials is defined as the difference 
between the actual state vector and the predicted state vector. 

The times represented on the x-axis of the two figures is in term of seconds since 
the beginning of the calendar year accordance with the time step of each state vector. 
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The Clohessy-Wiltshire propagator predicted sinusoidal variations as seen in 
Figures 5.3 (a) and (b). The cause of the sinusoidal motion is a result of the initial 
assumption in deriving the Clohessy-Wiltshire Equations that the target orbit is circular. 
In reality the target orbit is more nearly elliptical, thereby introducing an error. This fact 
was discussed in Chapter III, Section B. 

Additionally, large errors caused by thruster firings can be seen in both Figures 
5.3 and 5.4. The largest of these occur at identical times with the Encke, Cowell and 
Clohessy-Wiltshire propagators, but notice that the Clohessy-Wiltshire is more 
susceptible to small firings judging form the relatively large number of spikes occurring 
in Figure 5.3 compared to those in Figure 5.2. 



-X Coordinate 
-Y Coordinate 
~Z Coordinate 


Figure 5.3 (a) Position Differential of CW Propagator for a 3.84 Time Step 
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Figure 5.3 (b) Velocity Differential of CW Propagator for 3.84 Time Step 

It is easier to view the separate differences in the three propagators when they are 
compared with respect to Position and Velocity Components. This type of comparison 
also makes it possible to see the differences produced by adding perturbing accelerations 
to the Clohessy-Wiltshire propagator. For the purpose of this analysis the gravity model 
is applied at two different resolutions to determine whether or not the number of tesserals 
and sectorials is significant. Also the atmospheric model and gravity models are applied 
both separately and together. 

Additionally all propagators demonstrated a sinusoidal nature in the differentials. 
The commonality tends to suggest that it is the data set that is the cause vise the 
propagators. Possibly the error is a result of the manner in which the measurements were 
taken. An accelerometer, for example, can demonstrate a slight sinusoidal motion not 
caused by external accelerations. 







Figure 5.4 displays the X Coordinate Differential. This graph denotes the 
difference between the predicted state vector and actual state vector for the same time 
verses increasing propagation time steps. Each differential is the average result of 
propagating all 3238 state vectors the appropriate time step and taking the difference. As 
expected, the Cowell propagator shows the smallest differential even at the largest time 
step. The Encke propagator parallels the Cowell’s performance with error most likely 
due to its assumption of perfect elliptical motion. The Clohessy-Wiltshire shows the 
worst performance. The atmospheric model produced no apparent effect which really is 
not that surprising considering the low density of the atmosphere coupled with the short 
time steps consider. The fact that the addition of the gravitational model resulted in a 
worse solution not a better lends credence to the possibility that the data set contains 
processed state vectors not raw state vectors. 

In Figure 5.5 the Y Coordinate Differential is shown. All propagators show 
almost identical error, including the end point of the unperturbed Clohessy-Wiltshire. 
This may appear unusual until the magnitude of the Delta’s is examined, which at their 
worst never exceed 0.3 meters. The small magnitude of this error is because the motion 
of both the Orbiter and Target spacecraft are in very nearly the same orbital plane. 
Considering the small value of the input component the error is proportionately small. 

The Z Coordinate Differential, shown in Figure 5.6, shows predictable results. 
Again the Cowell and Encke propagators are essentially error free. The addition of 
gravitational and atmospheric drag perturbations in the Clohessy-Wiltshire produce a 
negligible effect. The reason for the negligible difference when the perturbing forces are 
added is most likely due to the magnitude of the velocity in the z-direction. Since 
atmospheric drag in a specific direction is a function of the velocity of the spacecraft in 
that direction, and because the x-axis of the coordinate system used is defined by the 
spacecraft’s velocity vector, drag can be expected to have a greater effect in the x- 
direction than in the z. 
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igure 5.4 X Coordinate Differential 
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Figure 5.5 Y Coordinate Differential 
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Over the given time interval, the velocity components never exceed 7.5 km/s 
while the velocity differentials are all less than 2.5 m/s. Figures 5.7, 5.8 and 5.9 show the 
Velocity Differentials produced. In all three cases the Cowell propagator displays the 
least error. The Encke propagator is generally the second most accurate, except in the 
case of the Y Component where the Encke incurs it’s greatest relative error. The source 
of this error can be traced back to the initial assumption made in deriving the “f and g” 
functions, which is that all motion is assumed to occur in the orbital plane, also known as 
the x and z plane. Therefore the variation in the component is a measure of the changes 
in out of plane relative velocity of the Orbiter and Target occurring over the time interval. 
This effect probably would be seen in the Y Position Component if the time interval had 
been greater. 

E. COMPUTATION TIME 

The computation times required for the each propagator can be seen in Table 5.1. 
They were determined by reading the computer’s, biostimer at program start and finish 
and dividing the difference by the number of vectors propagated. As expected, the 
Clohessy-Wiltshire proved the quickest with Encke next and the Cowell taking the 
longest time. Adding the perturbation models significantly increase the computation time 
of the Clohessy-Wiltshire, with little, if any, increase in accuracy. In fact, Encke 
propagator without any perturbations proved to be more accurate with the second shortest 
computation time. 
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Figure 5.8 Y Velocity Differentia] 







Figure 5.9 Z Velocity Differential 
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Table 5.1 Propagator Computation Times 
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VI. CONCLUSIONS 


The analysis performed by this thesis has attempted to determine which of three 
representative orbital propagators demonstrates the best performance. This performance 
determination is based not only on the accuracy of the propagator, but also on the 
computation time the propagator requires. Coupled with both of these criteria is the 
question of whether a propagator should apply the extensive models required for 
perturbing accelerations. Is the gain in accuracy worth the increased time required to 
compute the results. 

A. THE RESULTS 

As expected the Cowell propagator demonstrated the most accurate result while 
the Clohessy-Wiltshire required the least computational effort. As feared the perturbing 
accelerations produced a less accurate solution. This should be impossible, therefore the 
data set used must have been of processed data. This fact does not completely discount 
the project findings because a relative performance between the propagators can be 
inferred based upon the unperturbed results. 

The greatest performance appears to be demonstrated by the Encke propagator 
which without the addition of perturbing accelerations is essentially the “f and g” 
functions. The Encke propagator demonstrated near Cowell accuracy with near 
Clohessy-Wiltshire computation times. Only in propagation of the y velocity component 
did the Clohessy-Wiltshire demonstrate superior performance to that of the Encke 
methodology. 
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B. ADDITIONAL WORK 


In order to fully complete this project a number of things need to be 
accomplished. 

1. Raw Data 

The entire thesis should be repeated using the raw state vectors from a single 
sensing source. This sensing source need not be specified, but the data set used must not 
have been previously subjected to a propagator in order to smooth the data set as was 
probably done in this case. The C++ code developed in this thesis is already set up to 
perform this analysis whenever an appropriate data set can be found. 

2. Encke Propagator 

While the Encke Propagator worked extremely well without perturbing 
accelerations, this project still requires a working Encke propagator with perturbing 
accelerations in order to determine the potential benefits of adding these. The core work 
and theory exists, but the difficulties in compiling and executing the C++ code must be 
resolved. It is not anticipated that the addition of perturbing accelerations will result in 
any significant improvement in the algorithms performance for the short time interval 
calculations. This is especially true if one considers the fact that addition of the 
perturbing accelerations significantly increases computation time. 

3. Repeated Iterations 

The propagators should also be modified to check for the accuracy of the 
algorithms over longer propagation time in order to define the propagator limits. This 
could probably be done by breaking the propagation time into a series of sufficiently 
small time steps to maximize performance. This is identical to the method used in the 
Cowell propagator. In this method the state vector would be propagated over a series of 
time steps with the resulting state vector being used in successive iterations. The stability 
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of the algorithms as well as the additive effects of error should be critically examined. 
Should an operational algorithm be developed using one of the propagators reviewed in 
this thesis, such an analysis should be done in the selection process. 

4. Different Orbits 

This and the above projects should also be repeated for different orbits to 
determine which performs best under those circumstances. More detailed knowledge is 
needed of the limits at which the perturbing forces are required. In this thesis, some of 
the possible perturbations were discarded purely based on logical assumptions. These 
could be examined. Additionally the perturbations used probably have limits which 
should be examined. 
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